/*
 * DiagonalIncompleteLUFactorization.cs
 * 
 * Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

using System;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;

namespace dnAnalytics.LinearAlgebra.Solvers.Preconditioners
{
    /// <summary>
    /// Creates a diagonal incomplete LU preconditioner.
    /// </summary>
    /// <remarks>
    /// The D-Incomplete LU algorithm was taken from: <br/>
    /// Templates for the solution of linear systems: Building blocks
    /// for iterative methods <br/>
    /// Richard Barrett, Michael Berry, Tony F. Chan, James Demmel,
    /// June M. Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo,
    /// Charles Romine and Henk van der Vorst 
    /// <br/>
    /// Url: <a href="http://www.netlib.org/templates/Templates.html">http://www.netlib.org/templates/Templates.html</a>
    /// <br/>
    /// Algorithm is described in Chapter 3, section 3.4, page 43
    /// </remarks>
    /// <include file='../../../../../examples/trunk/LinearAlgebra/Solvers/Preconditioners/DiagonalIluFactorization.xml' path='example'/> 
    public sealed class DiagonalIncompleteLUFactorization : IPreconditioner
    {
        /// <summary>
        /// The original coefficient matrix.
        /// </summary>
        private SparseMatrix mCoefficientMatrix;
        /// <summary>
        /// The decomposed upper triangular matrix.
        /// </summary>
        private SparseMatrix mL;
        /// <summary>
        /// The decomposed lower triangular matrix.
        /// </summary>
        private SparseMatrix mU;
        /// <summary>
        /// The decomposed matrix diagonal values.
        /// </summary>
        private double[] mD;
        /// <summary>
        /// A flag used to indicate if the preconditioner has been created, i.e.
        /// if the <b>mLuDecomposition</b> matrix has been filled with the correct
        /// values.
        /// </summary>
        private bool mWasPreconditionerCreated;

        /// <summary>
        /// Creates a new <b>DiagonalIncompleteLUFactorization</b> instance.
        /// </summary>
        /// <param name="coefficientMatrix">The <see cref="SparseMatrix"/> upon which this
        /// preconditioner is based.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="coefficientMatrix"/> is <b>null</b>.</exception>
        /// <exception cref="NotSquareMatrixException">If <paramref name="coefficientMatrix"/> is not a square matrix.</exception>
        public DiagonalIncompleteLUFactorization(Matrix coefficientMatrix)
        {
            if (coefficientMatrix == null)
            {
                throw new ArgumentNullException("coefficientMatrix", Strings.NullParameterException);
            }

            if (coefficientMatrix.Rows != coefficientMatrix.Columns)
            {
                throw new NotSquareMatrixException();
            }

            mCoefficientMatrix = new SparseMatrix(coefficientMatrix);
        }

        /// <summary>
        /// Returns the decomposed upper triangular matrix.
        /// </summary>
        /// <returns>The upper triangular matrix.</returns>
        internal SparseMatrix UpperTriagonalMatrix()
        {
            return mU;
        }

        /// <summary>
        /// Returns the decomposed lower triangular matrix.
        /// </summary>
        /// <returns>The lower triangular matrix.</returns>
        internal SparseMatrix LowerTriagonalMatrix()
        {
            return mL;
        }

        /// <summary>
        /// Returns the decomposed matrix diagonal.
        /// </summary>
        /// <returns>The matrix diagonal.</returns>
        internal SparseMatrix DiagonalMatrix()
        {
            SparseMatrix result = new SparseMatrix(mD.Length);
            for (int i = 0; i < mD.Length; i++)
            {
                result[i, i] = 1 / mD[i];
            }
            return result;
        }

        /// <summary>
        /// Initializes the preconditioner and loads the internal data structures.
        /// </summary>
        public void CreatePreconditioner()
        {
            /* 
            Create the preconditioning matrix. We need D, L and U
            (the Diagonal matrix, strictly Lower and strictly Upper triangular
            matrices). 
            The diagonal matrix is a single array of values created accoording to
            for (int i = 0; i < matrix.RowLength; i++)
            {
                d_ii = a_ii
            }
            for (int i = 0; i < matrix.RowLength; i++)
            {
                d_ii = 1/d_ii // NOTE THIS CAN NOT BE DONE EARLIER DUE TO THE NEXT STEP
                for (int j = i ; j < matrix.RowLength; j++)
                {
                    if ((a_ij != 0.0) && (a_ji != 0.0))
                    {
                        d_jj = d_jj - a_ji * d_ii * a_ij
                    }
                }
            }
            */

            if (!mWasPreconditionerCreated)
            {

                mL = new SparseMatrix(mCoefficientMatrix.Rows);
                mU = new SparseMatrix(mCoefficientMatrix.Rows);
                for (int i = 0; i < mCoefficientMatrix.Rows; i++)
                {
                    for (int j = 0; j < mCoefficientMatrix.Columns; j++)
                    {
                        // Copy the upper triagonal matrix to mU
                        if (i < j)
                        {
                            mU[i, j] = mCoefficientMatrix[i, j];
                        }
                        // Copy the lower triagonal matrix to mL
                        if (i > j)
                        {
                            mL[i, j] = mCoefficientMatrix[i, j];
                        }
                    }
                }

                // Create the diagonal matrix
                mD = new double[mCoefficientMatrix.Rows];
                for (int i = 0; i < mCoefficientMatrix.Rows; i++)
                {
                    mD[i] = mCoefficientMatrix[i, i];
                }

                for (int i = 0; i < mCoefficientMatrix.Rows; i++)
                {
                    mD[i] = 1 / mD[i];
                    for (int j = i + 1; j < mCoefficientMatrix.Rows; j++)
                    {
                        if ((System.Math.Abs(mCoefficientMatrix[i, j]) > Double.Epsilon) &&
                            (System.Math.Abs(mCoefficientMatrix[j, i]) > Double.Epsilon))
                        {
                            mD[j] = mD[j] - mCoefficientMatrix[j, i] * mD[i] * mCoefficientMatrix[i, j];
                        }
                    }
                }
                mWasPreconditionerCreated = true;
            }
        }

        /// <summary>
        /// Solves the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="vector">The right hand side vector.</param>
        /// <param name="result">The left hand side vector.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="vector"/> is <b>null</b>.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If <paramref name="vector"/> and <paramref name="result"/> do not
        /// have the same size.</para>
        /// <para> - or -</para>
        /// <para>If the size of <paramref name="result"/> is different the number of rows of the coefficient matrix.</para></exception>
        public void Solve(Vector vector, Vector result)
        {
            /* 
            Allow to solve system Mx = y
            via:
            (D + L)z = y
            (I + D^-1 U) x = z
            Which gives
            for (int i = 1; i < matrix.RowLength; i++)
            {
                z_i = l_ii^-1 * (y_i - SUM_(j<i) l_ij * z_j)
            }
            for (int i = matrix.RowLength - 1; i > -1; i--)
            {
                x_i = u_ii^-1 * (z_i - SUM_(j > i) u_ij * x_j)
            }
            */

            if (vector == null)
            {
                throw new ArgumentNullException("vector", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Count != vector.Count) || (result.Count != mCoefficientMatrix.Rows))
            {
                throw new NotConformableException(Strings.ParametersNotConformable);
            }

            if (!mWasPreconditionerCreated)
            {
                CreatePreconditioner();
            }

            /*
            Solve:
            (D + L)z = y
            Which gives
            for (int i = 1; i < matrix.RowLength; i++)
            {
                z_i = d_ii^-1 * (y_i - SUM_(j<i) l_ij * z_j)
            }
            */
            Vector rowVector = VectorBuilder.CreateVector(mCoefficientMatrix.Rows, VectorType.Sparse);
            for (int i = 0; i < mCoefficientMatrix.Rows; i++)
            {
                mL.Row(i, rowVector);

                double sum = rowVector.DotProduct(result);
                /*for (int j = 0; j < i; j++)
                {
                    sum += rowVector[j] * result[j];
                }*/

                result[i] = mD[i] * (vector[i] - sum);
            }

            /*    
            Solve:
            (D + U) x = z
            Which gives
            for (int i = matrix.RowLength - 1; i > -1; i--)
            {
                x_i = z_i - d_ii^-1 * SUM_(j > i) u_ij * x_j
            }
            */
            for (int i = mCoefficientMatrix.Rows - 1; i > -1; i--)
            {
                mU.Row(i, rowVector);

                double sum = rowVector.DotProduct(result);
                /*for (int j = mCoefficientMatrix.Rows - 1; j > i; j--)
                {
                    sum += rowVector[j] * result[j];
                }*/

                result[i] = result[i] - mD[i] * sum;
            }
        }

        /// <summary>
        /// Solves the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <remarks>
        /// When using the <b>Solve</b> method consider carefully if the
        /// left hand side vector can be reused. If so it would be better to 
        /// use the overloaded <b>Solve</b> method.
        /// </remarks>
        /// <param name="vector">The right hand side vector.</param>
        /// <returns>The left hand side vector.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="vector"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If the size of <paramref name="vector"/> is different the number of rows of the coefficient matrix.</para></exception>
        public Vector Solve(Vector vector)
        {
            if (vector == null)
            {
                throw new ArgumentNullException("vector", Strings.NullParameterException);
            }

            if (vector.Count != mCoefficientMatrix.Columns)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            Vector result = VectorBuilder.CreateVector(vector.Count, VectorType.Dense);
            Solve(vector, result);
            return result;
        }
    }
}
